home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 26 / Cream of the Crop 26.iso / os2 / octa209s.zip / octave-2.09 / libcruft / lapack / ztrsyl.f < prev    next >
Text File  |  1996-07-19  |  12KB  |  369 lines
  1.       SUBROUTINE ZTRSYL( TRANA, TRANB, ISGN, M, N, A, LDA, B, LDB, C,
  2.      $                   LDC, SCALE, INFO )
  3. *
  4. *  -- LAPACK routine (version 2.0) --
  5. *     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
  6. *     Courant Institute, Argonne National Lab, and Rice University
  7. *     March 31, 1993
  8. *
  9. *     .. Scalar Arguments ..
  10.       CHARACTER          TRANA, TRANB
  11.       INTEGER            INFO, ISGN, LDA, LDB, LDC, M, N
  12.       DOUBLE PRECISION   SCALE
  13. *     ..
  14. *     .. Array Arguments ..
  15.       COMPLEX*16         A( LDA, * ), B( LDB, * ), C( LDC, * )
  16. *     ..
  17. *
  18. *  Purpose
  19. *  =======
  20. *
  21. *  ZTRSYL solves the complex Sylvester matrix equation:
  22. *
  23. *     op(A)*X + X*op(B) = scale*C or
  24. *     op(A)*X - X*op(B) = scale*C,
  25. *
  26. *  where op(A) = A or A**H, and A and B are both upper triangular. A is
  27. *  M-by-M and B is N-by-N; the right hand side C and the solution X are
  28. *  M-by-N; and scale is an output scale factor, set <= 1 to avoid
  29. *  overflow in X.
  30. *
  31. *  Arguments
  32. *  =========
  33. *
  34. *  TRANA   (input) CHARACTER*1
  35. *          Specifies the option op(A):
  36. *          = 'N': op(A) = A    (No transpose)
  37. *          = 'C': op(A) = A**H (Conjugate transpose)
  38. *
  39. *  TRANB   (input) CHARACTER*1
  40. *          Specifies the option op(B):
  41. *          = 'N': op(B) = B    (No transpose)
  42. *          = 'C': op(B) = B**H (Conjugate transpose)
  43. *
  44. *  ISGN    (input) INTEGER
  45. *          Specifies the sign in the equation:
  46. *          = +1: solve op(A)*X + X*op(B) = scale*C
  47. *          = -1: solve op(A)*X - X*op(B) = scale*C
  48. *
  49. *  M       (input) INTEGER
  50. *          The order of the matrix A, and the number of rows in the
  51. *          matrices X and C. M >= 0.
  52. *
  53. *  N       (input) INTEGER
  54. *          The order of the matrix B, and the number of columns in the
  55. *          matrices X and C. N >= 0.
  56. *
  57. *  A       (input) COMPLEX*16 array, dimension (LDA,M)
  58. *          The upper triangular matrix A.
  59. *
  60. *  LDA     (input) INTEGER
  61. *          The leading dimension of the array A. LDA >= max(1,M).
  62. *
  63. *  B       (input) COMPLEX*16 array, dimension (LDB,N)
  64. *          The upper triangular matrix B.
  65. *
  66. *  LDB     (input) INTEGER
  67. *          The leading dimension of the array B. LDB >= max(1,N).
  68. *
  69. *  C       (input/output) COMPLEX*16 array, dimension (LDC,N)
  70. *          On entry, the M-by-N right hand side matrix C.
  71. *          On exit, C is overwritten by the solution matrix X.
  72. *
  73. *  LDC     (input) INTEGER
  74. *          The leading dimension of the array C. LDC >= max(1,M)
  75. *
  76. *  SCALE   (output) DOUBLE PRECISION
  77. *          The scale factor, scale, set <= 1 to avoid overflow in X.
  78. *
  79. *  INFO    (output) INTEGER
  80. *          = 0: successful exit
  81. *          < 0: if INFO = -i, the i-th argument had an illegal value
  82. *          = 1: A and B have common or very close eigenvalues; perturbed
  83. *               values were used to solve the equation (but the matrices
  84. *               A and B are unchanged).
  85. *
  86. *  =====================================================================
  87. *
  88. *     .. Parameters ..
  89.       DOUBLE PRECISION   ONE
  90.       PARAMETER          ( ONE = 1.0D+0 )
  91. *     ..
  92. *     .. Local Scalars ..
  93.       LOGICAL            NOTRNA, NOTRNB
  94.       INTEGER            J, K, L
  95.       DOUBLE PRECISION   BIGNUM, DA11, DB, EPS, SCALOC, SGN, SMIN,
  96.      $                   SMLNUM
  97.       COMPLEX*16         A11, SUML, SUMR, VEC, X11
  98. *     ..
  99. *     .. Local Arrays ..
  100.       DOUBLE PRECISION   DUM( 1 )
  101. *     ..
  102. *     .. External Functions ..
  103.       LOGICAL            LSAME
  104.       DOUBLE PRECISION   DLAMCH, ZLANGE
  105.       COMPLEX*16         ZDOTC, ZDOTU
  106.       EXTERNAL           LSAME, DLAMCH, ZLANGE, ZDOTC, ZDOTU
  107. *     ..
  108. *     .. External Subroutines ..
  109.       EXTERNAL           DLABAD, XERBLA, ZDSCAL
  110. *     ..
  111. *     .. Intrinsic Functions ..
  112.       INTRINSIC          ABS, DBLE, DCMPLX, DCONJG, DIMAG, MAX, MIN
  113. *     ..
  114. *     .. Executable Statements ..
  115. *
  116. *     Decode and Test input parameters
  117. *
  118.       NOTRNA = LSAME( TRANA, 'N' )
  119.       NOTRNB = LSAME( TRANB, 'N' )
  120. *
  121.       INFO = 0
  122.       IF( .NOT.NOTRNA .AND. .NOT.LSAME( TRANA, 'T' ) .AND. .NOT.
  123.      $    LSAME( TRANA, 'C' ) ) THEN
  124.          INFO = -1
  125.       ELSE IF( .NOT.NOTRNB .AND. .NOT.LSAME( TRANB, 'T' ) .AND. .NOT.
  126.      $         LSAME( TRANB, 'C' ) ) THEN
  127.          INFO = -2
  128.       ELSE IF( ISGN.NE.1 .AND. ISGN.NE.-1 ) THEN
  129.          INFO = -3
  130.       ELSE IF( M.LT.0 ) THEN
  131.          INFO = -4
  132.       ELSE IF( N.LT.0 ) THEN
  133.          INFO = -5
  134.       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
  135.          INFO = -7
  136.       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
  137.          INFO = -9
  138.       ELSE IF( LDC.LT.MAX( 1, M ) ) THEN
  139.          INFO = -11
  140.       END IF
  141.       IF( INFO.NE.0 ) THEN
  142.          CALL XERBLA( 'ZTRSYL', -INFO )
  143.          RETURN
  144.       END IF
  145. *
  146. *     Quick return if possible
  147. *
  148.       IF( M.EQ.0 .OR. N.EQ.0 )
  149.      $   RETURN
  150. *
  151. *     Set constants to control overflow
  152. *
  153.       EPS = DLAMCH( 'P' )
  154.       SMLNUM = DLAMCH( 'S' )
  155.       BIGNUM = ONE / SMLNUM
  156.       CALL DLABAD( SMLNUM, BIGNUM )
  157.       SMLNUM = SMLNUM*DBLE( M*N ) / EPS
  158.       BIGNUM = ONE / SMLNUM
  159.       SMIN = MAX( SMLNUM, EPS*ZLANGE( 'M', M, M, A, LDA, DUM ),
  160.      $       EPS*ZLANGE( 'M', N, N, B, LDB, DUM ) )
  161.       SCALE = ONE
  162.       SGN = ISGN
  163. *
  164.       IF( NOTRNA .AND. NOTRNB ) THEN
  165. *
  166. *        Solve    A*X + ISGN*X*B = scale*C.
  167. *
  168. *        The (K,L)th block of X is determined starting from
  169. *        bottom-left corner column by column by
  170. *
  171. *            A(K,K)*X(K,L) + ISGN*X(K,L)*B(L,L) = C(K,L) - R(K,L)
  172. *
  173. *        Where
  174. *                    M                        L-1
  175. *          R(K,L) = SUM [A(K,I)*X(I,L)] +ISGN*SUM [X(K,J)*B(J,L)].
  176. *                  I=K+1                      J=1
  177. *
  178.          DO 30 L = 1, N
  179.             DO 20 K = M, 1, -1
  180. *
  181.                SUML = ZDOTU( M-K, A( K, MIN( K+1, M ) ), LDA,
  182.      $                C( MIN( K+1, M ), L ), 1 )
  183.                SUMR = ZDOTU( L-1, C( K, 1 ), LDC, B( 1, L ), 1 )
  184.                VEC = C( K, L ) - ( SUML+SGN*SUMR )
  185. *
  186.                SCALOC = ONE
  187.                A11 = A( K, K ) + SGN*B( L, L )
  188.                DA11 = ABS( DBLE( A11 ) ) + ABS( DIMAG( A11 ) )
  189.                IF( DA11.LE.SMIN ) THEN
  190.                   A11 = SMIN
  191.                   DA11 = SMIN
  192.                   INFO = 1
  193.                END IF
  194.                DB = ABS( DBLE( VEC ) ) + ABS( DIMAG( VEC ) )
  195.                IF( DA11.LT.ONE .AND. DB.GT.ONE ) THEN
  196.                   IF( DB.GT.BIGNUM*DA11 )
  197.      $               SCALOC = ONE / DB
  198.                END IF
  199.                X11 = ( VEC*DCMPLX( SCALOC ) ) / A11
  200. *
  201.                IF( SCALOC.NE.ONE ) THEN
  202.                   DO 10 J = 1, N
  203.                      CALL ZDSCAL( M, SCALOC, C( 1, J ), 1 )
  204.    10             CONTINUE
  205.                   SCALE = SCALE*SCALOC
  206.                END IF
  207.                C( K, L ) = X11
  208. *
  209.    20       CONTINUE
  210.    30    CONTINUE
  211. *
  212.       ELSE IF( .NOT.NOTRNA .AND. NOTRNB ) THEN
  213. *
  214. *        Solve    A' *X + ISGN*X*B = scale*C.
  215. *
  216. *        The (K,L)th block of X is determined starting from
  217. *        upper-left corner column by column by
  218. *
  219. *            A'(K,K)*X(K,L) + ISGN*X(K,L)*B(L,L) = C(K,L) - R(K,L)
  220. *
  221. *        Where
  222. *                   K-1                         L-1
  223. *          R(K,L) = SUM [A'(I,K)*X(I,L)] + ISGN*SUM [X(K,J)*B(J,L)]
  224. *                   I=1                         J=1
  225. *
  226.          DO 60 L = 1, N
  227.             DO 50 K = 1, M
  228. *
  229.                SUML = ZDOTC( K-1, A( 1, K ), 1, C( 1, L ), 1 )
  230.                SUMR = ZDOTU( L-1, C( K, 1 ), LDC, B( 1, L ), 1 )
  231.                VEC = C( K, L ) - ( SUML+SGN*SUMR )
  232. *
  233.                SCALOC = ONE
  234.                A11 = DCONJG( A( K, K ) ) + SGN*B( L, L )
  235.                DA11 = ABS( DBLE( A11 ) ) + ABS( DIMAG( A11 ) )
  236.                IF( DA11.LE.SMIN ) THEN
  237.                   A11 = SMIN
  238.                   DA11 = SMIN
  239.                   INFO = 1
  240.                END IF
  241.                DB = ABS( DBLE( VEC ) ) + ABS( DIMAG( VEC ) )
  242.                IF( DA11.LT.ONE .AND. DB.GT.ONE ) THEN
  243.                   IF( DB.GT.BIGNUM*DA11 )
  244.      $               SCALOC = ONE / DB
  245.                END IF
  246. *
  247.                X11 = ( VEC*DCMPLX( SCALOC ) ) / A11
  248. *
  249.                IF( SCALOC.NE.ONE ) THEN
  250.                   DO 40 J = 1, N
  251.                      CALL ZDSCAL( M, SCALOC, C( 1, J ), 1 )
  252.    40             CONTINUE
  253.                   SCALE = SCALE*SCALOC
  254.                END IF
  255.                C( K, L ) = X11
  256. *
  257.    50       CONTINUE
  258.    60    CONTINUE
  259. *
  260.       ELSE IF( .NOT.NOTRNA .AND. .NOT.NOTRNB ) THEN
  261. *
  262. *        Solve    A'*X + ISGN*X*B' = C.
  263. *
  264. *        The (K,L)th block of X is determined starting from
  265. *        upper-right corner column by column by
  266. *
  267. *            A'(K,K)*X(K,L) + ISGN*X(K,L)*B'(L,L) = C(K,L) - R(K,L)
  268. *
  269. *        Where
  270. *                    K-1
  271. *           R(K,L) = SUM [A'(I,K)*X(I,L)] +
  272. *                    I=1
  273. *                           N
  274. *                     ISGN*SUM [X(K,J)*B'(L,J)].
  275. *                          J=L+1
  276. *
  277.          DO 90 L = N, 1, -1
  278.             DO 80 K = 1, M
  279. *
  280.                SUML = ZDOTC( K-1, A( 1, K ), 1, C( 1, L ), 1 )
  281.                SUMR = ZDOTC( N-L, C( K, MIN( L+1, N ) ), LDC,
  282.      $                B( L, MIN( L+1, N ) ), LDB )
  283.                VEC = C( K, L ) - ( SUML+SGN*DCONJG( SUMR ) )
  284. *
  285.                SCALOC = ONE
  286.                A11 = DCONJG( A( K, K )+SGN*B( L, L ) )
  287.                DA11 = ABS( DBLE( A11 ) ) + ABS( DIMAG( A11 ) )
  288.                IF( DA11.LE.SMIN ) THEN
  289.                   A11 = SMIN
  290.                   DA11 = SMIN
  291.                   INFO = 1
  292.                END IF
  293.                DB = ABS( DBLE( VEC ) ) + ABS( DIMAG( VEC ) )
  294.                IF( DA11.LT.ONE .AND. DB.GT.ONE ) THEN
  295.                   IF( DB.GT.BIGNUM*DA11 )
  296.      $               SCALOC = ONE / DB
  297.                END IF
  298. *
  299.                X11 = ( VEC*DCMPLX( SCALOC ) ) / A11
  300. *
  301.                IF( SCALOC.NE.ONE ) THEN
  302.                   DO 70 J = 1, N
  303.                      CALL ZDSCAL( M, SCALOC, C( 1, J ), 1 )
  304.    70             CONTINUE
  305.                   SCALE = SCALE*SCALOC
  306.                END IF
  307.                C( K, L ) = X11
  308. *
  309.    80       CONTINUE
  310.    90    CONTINUE
  311. *
  312.       ELSE IF( NOTRNA .AND. .NOT.NOTRNB ) THEN
  313. *
  314. *        Solve    A*X + ISGN*X*B' = C.
  315. *
  316. *        The (K,L)th block of X is determined starting from
  317. *        bottom-left corner column by column by
  318. *
  319. *           A(K,K)*X(K,L) + ISGN*X(K,L)*B'(L,L) = C(K,L) - R(K,L)
  320. *
  321. *        Where
  322. *                    M                          N
  323. *          R(K,L) = SUM [A(K,I)*X(I,L)] + ISGN*SUM [X(K,J)*B'(L,J)]
  324. *                  I=K+1                      J=L+1
  325. *
  326.          DO 120 L = N, 1, -1
  327.             DO 110 K = M, 1, -1
  328. *
  329.                SUML = ZDOTU( M-K, A( K, MIN( K+1, M ) ), LDA,
  330.      $                C( MIN( K+1, M ), L ), 1 )
  331.                SUMR = ZDOTC( N-L, C( K, MIN( L+1, N ) ), LDC,
  332.      $                B( L, MIN( L+1, N ) ), LDB )
  333.                VEC = C( K, L ) - ( SUML+SGN*DCONJG( SUMR ) )
  334. *
  335.                SCALOC = ONE
  336.                A11 = A( K, K ) + SGN*DCONJG( B( L, L ) )
  337.                DA11 = ABS( DBLE( A11 ) ) + ABS( DIMAG( A11 ) )
  338.                IF( DA11.LE.SMIN ) THEN
  339.                   A11 = SMIN
  340.                   DA11 = SMIN
  341.                   INFO = 1
  342.                END IF
  343.                DB = ABS( DBLE( VEC ) ) + ABS( DIMAG( VEC ) )
  344.                IF( DA11.LT.ONE .AND. DB.GT.ONE ) THEN
  345.                   IF( DB.GT.BIGNUM*DA11 )
  346.      $               SCALOC = ONE / DB
  347.                END IF
  348. *
  349.                X11 = ( VEC*DCMPLX( SCALOC ) ) / A11
  350. *
  351.                IF( SCALOC.NE.ONE ) THEN
  352.                   DO 100 J = 1, N
  353.                      CALL ZDSCAL( M, SCALOC, C( 1, J ), 1 )
  354.   100             CONTINUE
  355.                   SCALE = SCALE*SCALOC
  356.                END IF
  357.                C( K, L ) = X11
  358. *
  359.   110       CONTINUE
  360.   120    CONTINUE
  361. *
  362.       END IF
  363. *
  364.       RETURN
  365. *
  366. *     End of ZTRSYL
  367. *
  368.       END
  369.